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To the present day, the Beenakker-Mazur (BM) method is the most comprehensive statistical 
physics approach to the calculation of short-time transport properties of colloidal suspensions. A 
revised version of the BM method with an improved treatment of hydrodynamic interactions is pre¬ 
sented and evaluated regarding the rotational short-time self-diffusion coefficient, D r , of suspensions 
of charged particles interacting by a hard-sphere plus screened Coulomb (Yukawa) pair potential. 

To assess the accuracy of the method, elaborate simulations of D r have been performed, covering 
a broad range of interaction parameters and particle concentrations. The revised BM method is 
compared in addition with results by a simplifying pairwise additivity (PA) method in which the 
hydrodynamic interactions are treated on a two-body level. The static pair correlation functions re¬ 
quired as input to both theoretical methods are calculated using the Rogers-Young integral equation 
scheme. While the revised BM method reproduces the general trends of the simulation results, it 
systematically and significantly underestimates the rotational diffusion coefficient. The PA method 
agrees well with the simulation data at lower volume fractions, but at higher concentrations D r is 
likewise underestimated. For a fixed value of the pair potential at mean particle distance comparable 
to the thermal energy, D r increases strongly with increasing Yukawa potential screening parameter. 
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I. INTRODUCTION 

Short-time transport properties of colloidal suspensions such as translational self- and collective diffusion coefficients, 
hydrodynamic function and high-frequency viscosity have been the subject of numerous experimental studies [ 1 — 
7]. These studies have been accompanied over the past years by computer simulation studies (see, e.g. [ 8 , 9]) 
and theoretical works [10-19]. While short-time transport properties are expressible theoretically as rather simple 
equilibrium averages invoking hydrodynamic mobility tensors, the difficulty in their actual calculation arises from the 
long-ranged many-body hydrodynamic interactions (His) between the particles. The slowing influence of the His is 
particularly pronounced when particles are in relative motion close to each other. 

There exist two major solution schemes which have been used for the calculation of short-time properties. The 
first one is the so-called pairwise additivity (PA) approximation. In its most complete version, all two-body His 
contributions are accounted for including lubrication terms, but three-body and higher order interaction contributions 
are disregarded. By construction, its application range is usually limited to semi-dilute systems, and here in particular 
to charge-stabilized suspensions where near-contact configurations of particles are statistically unlikely [9] . Regarding 
short-time self-diffusion coefficients and the high-frequency viscosity, however, the PA method can be profitably used 
also at higher concentrations, owing to the steep decay, with increasing inter-particle distances, of the hydrodynamic 
mobility tensors associated with these quantities. 

Different from the PA scheme, the semi-analytical method of calculating short-time diffusion and viscosity properties 
by Beenakker-Mazur [13-1 ], commonly referred to as the £7 method, is applicable also to concentrated suspensions. 
It is a mean-field-type approximation method which accounts for many-body His contributions in the form of so- 
called ring self-correlation diagrams, without an account of lubrication effects. In its standard second-order version, 
the only required external input is the static structure factor S(q) as function of the scattering wavenumber q. A 
major short-coming of the original 5 7 method is its poor treatment of the translational short-time self-diffusion 
coefficient D t . This deficiency can be overcome by using a more accurate method for the self-part, e.g. by using the 
PA approximation result for D t at lower concentrated systems. The self-part corrected S"f scheme has been applied 
both to suspensions of neutral and charge-stabilized particles [9, 20, 21], for the calculation of the hydrodynamic 
and collective diffusion functions, and the high-frequency viscosity. It was applied recently also to suspensions of 
hydrodynamically structured particles [ 22 ] , and with additional approximations also to a binary hard-sphere mixture 
[23]. The predictions for these systems are decently good, with inaccuracies revealed at all concentrations. These 
inaccuracies can be partially attributed to the approximate treatment of the His in the £7 method, and partially to 
the invoked mean-field approximation. In recent work by Makuch and Cichocki [24], the approximation steps in the 
original derivation of the £7 method have been reduced, in particular by accounting for a large number of hydrodynamic 
multipoles in the truncated multipolar matrices which are extrapolated to infinite order. The observation that the 
revised ^7 method by Makuch and Cichocki, with its improved hydrodynamic mobility tensor treatment, has not 
resulted in a systematic improvement of the hydrodynamic function and high-frequency viscosity predictions points 
to a fortuitous cancellation of errors in the various approximation steps of the original method by Beenakker and 
Mazur. 

Experimental studies of the short-time rotational self-diffusion coefficient, D r , in concentrated suspensions are 
based on techniques which can distinguish different particle orientations. These methods include depolarized dynamic 
light scattering on optically anisotropic particles [3], nuclear magnetic resonance [25], time-resolved phosphorescence 
anisotropy [26-28] and polarized fluorescence recovery after photobleaching [29]. Most experimental work on rotational 
diffusion has been for monodisperse colloidal systems of neutral hard spheres and charge-stabilized particles. In 
addition, binary mixtures of charge-stabilized particles have been studied where one component is very dilute [26-28] . 

Simulation work on rotational diffusion has dealt so far mainly with monodisperse systems of non-permeable [21, 
30, 31] and permeable hard spheres [32]. Charge-stabilized systems have been addressed in few simulation studies 
only [21, 30], focused on low-salinity systems with weak electrostatic screening. Therefore, a systematic simulation 
study of rotational diffusion in charge-stabilized systems is still on demand. 

Short-time rotational self-diffusion in hard-splrere suspensions was studied theoretically by various groups using 
truncated hydrodynamic cluster expansions up to quadratic order in the particle volume fraction (f) [3, 33-35]. The 
high-precision result by Cichocki et al. [36], 

% = l-0.631<(>^0.726<)> 2 + 0(</> 3 ), (1) 

includes a lubrication correction for the three-body His contributions. Here, Dq = ksT/(8n^a 3 ) is the single-particle 
rotational diffusion coefficient of a no-slip sphere of hydrodynamic radius a. As shown in [21, 32], Eq. (1) describes 
simulation and experimental data remarkably well for volume fractions up to the freezing transition value 4>f = 0.494, 
indicating that higher-order virial coefficients are small or mutually cancel out. 
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Regarding rotational diffusion, the Sj method has been applied so far to hard-sphere suspensions only, in the work 
by Treloar and Masters [37] where approximations along the line of those introduced by Beenakker and Mazur have 
been made. It was not applied so far to rotational diffusion in charge-stabilized suspensions. However, these systems 
have been analyzed in the weak electrostatic screening regime using a simplified PA approximation approach based 
on a truncated inverse distance expansion of the two-body rotational mobility tensors, by considering in addition the 
leading-order long-distance hydrodynamic three-body term [27, 38, 39]. For charged particles with strong long-distance 
repulsion, the remarkable scaling relation, 
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with a r ~ 1.3 has been obtained by this simplifying approach. The comparison with Lattice-Boltzmann [30] and 
accelerated Stokesian dynamics simulation results [21] has shown that it applies accurately up to <j> ss 0.3. Different 
from Eq. (1) valid for neutral hard spheres, the scaling relation in Eq. (2) is not a second-order virial expansion result. 
It originates basically from the concentration scaling of the radius, r m , of next-neighbor shells in low-salinity 

systems [ 21 ]. 

The present article includes a comprehensive theoretical and simulation analysis of short-time rotational diffusion 
in fluid-like, charged particles suspensions whose static pair interactions are modeled by the hard-sphere plus repulsive 
Yukawa (HSY) pair potential. A revised version of the S'y method for D r by Treloar and Masters is evaluated for 
a broad range of volume fractions and reduced pair potential strengths, and two screening parameters characteristic 
for the weak and strong screening regimes, respectively. The accuracy of the revised 5 7 method for D r is assessed in 
the comparison with high-accuracy simulation results which we have obtained using a multipole simulation method 
encoded in the HYDROMULTPOLE simulation package [36]. The revised S'y method and simulation results are 
compared in addition with our predictions by a simplifying pairwise additivity (PA) method in which the rotational 
hydrodynamic mobility tensors are treated on the two-body level, without an additional long-distance truncation as 
made in earlier work. The radial distribution function (RDF), g(r), and static structure factor of the HSY model 
constituting the static input to the theoretical methods are calculated using the accurate Rogers-Young (RY) integral 
equation scheme. The present work is the first systematic theory-simulation study of rotational self-diffusion in the 
HSY model. 

The article is organized as follows: In Sec. II, we give the essentials of short-time rotational self-diffusion. Sec. Ill 
includes the description of the HSY model with employed interaction parameters, and a discussion of the RY radial 
distribution functions used in the analytical-theoretical calculation of D r . The revised Beenakker-Mazur method of 
calculating D r is explained in Sec. IV. In Secs. V and VI, respectively, the employed simulation and PA methods are 
described. Our results are presented in Sec. VII, and our finalizing conclusions in Sec. VIII. 


II. SHORT-TIME ROTATIONAL SELF-DIFFUSION COEFFICIENT 

We consider rotational diffusion in a fluid-state suspension of monodisperse spherical Brownian particles immersed 
in a structureless Newtonian solvent of shear viscosity 77 , on a coarse-grained Brownian time scale exceeding the 
rotational and translational momentum relaxation times Tg ~ tJ), respectively, by several orders of magnitude [3, 40]. 
On this scale, particles and fluid move quasi-inertia-free, and the fluid-mediated His are acting quasi-instantaneously. 
The configurational evolution of the particles is then governed by the generalized Smoluchowski equation [ 11-43] for 
the V-particle probability density function p (Ri,..., Rjv, ui,..., ujv, t) of the sphere center positions Ri,..., Rjv 
and orientations ui,..., ujv at time t. The associated low-Reynolds-number incompressible fluid flow is described by 
the linear stationary Stokes equation [44]. The hydrodynamic ingredients to the generalized Smoluchowski equation 
derived from the Stokes equation are the translational-rotational mobility tensors quantifying the linear relations, 
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between the forces and torques, Fj and T,, respectively, acting on the colloidal spheres, and the resulting translational 
and rotational particle velocities Uj and fi,;. For the uniformly assumed no-slip hydrodynamic surface boundary 
condition, the mobility tensors are independent of the particle orientations. 

In depolarized dynamic light scattering [3, 45], the short-time rotational self-diffusion coefficient of Brownian spheres 
is determined from the initial decay of the measurable orientational self-correlation function, for correlation times t 
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within Tg <C t -C 1 /Dq where particle orientations and positions have changed by very small amounts only, on the 
characteristic length scale of the suspension. For a concentrated isotropic suspension, D r can be computed as the 
ensemble average of the trace of the rotational mobility tensor [46], 
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where (• ■ ■) is an equilibrium ensemble average, and where the thermodynamic limit TV —^ oo at fixed particle 
concentration has been taken. It should be noted that D r as given in Eq. (5) is, for non-zero concentrations, different 
from the initial slope of the mean-squared displacement of the particle orientation unit vector fh(i). 


III. STATIC CORRELATIONS OF HSY PARTICLES 


The pair potential in the hard-sphere plus repulsive Yukawa (HSY) model is given by 

exp{— n{r — a)} 
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where 7 > 0 is the coupling parameter of the Yukawa-type potential part, <r = 2a the hard-core diameter, r the 
center-to-center distance between two spheres, k B the Boltzmann constant, and T the absolute temperature. The 
range of the HSY potential is set by the inverse of the screening parameter k > 0. In the infinite screening limit 
k —» 00 , or likewise for 7 = 0 , the hard-sphere potential is recovered. In the opposite limit k —> 0 of zero-screening, a 
one-component plasma-like system is described. 

Dispersions which can be described by the HSY model range from charge-stabilized suspensions of rigid colloidal 
spheres [47] to globular protein solutions [48] and dusty plasmas [ 19]. The HSY potential form is in general a 
good approximation to the state-dependent effective pair-potential between charged colloidal spheres. The latter is 
obtained from integrating out the degrees of freedom of the microions and solvent molecules. In many experimentally 
encountered suspensions, the short-ranged van der Waals attraction neglected in the HSY model is of no relevance, 
either since the electrostatic repulsion between the particles is strong enough to prevent near-contact configurations [7] , 
or the solvent dielectric constant nearly matches that of the suspended particles [50-53] , or the particles are sterically 
stabilized by grafted polymers [54], The complicated dependencies of the state-dependent potential parameters 7 
and k in charge-stabilized suspensions on the salt ion concentration, colloidal volume fraction, and bare and effective 
colloidal surface charges is the topic of on-going research that covers experiments, theory and computer simulations 
[55-66]. The present work is not concerned with a first-principles determination of the state-dependence of 7 and n 
which is also influenced by the specific surface electrochemistry of the dispersed particles. Instead, 7 and k are treated 
quite generally as individually variable parameters. Note further that the direct interactions in the HSY model are 
treated as pairwise additive. Non-pairwise additivity effects in the direct interaction of charged colloidal particles are 
usually quite small [67]. 

According to Eq. (6), the thermodynamic state of the HSY model, and likewise the RDF as a function of r/cr, are 
fully characterized by three independent non-dimensional parameters which can be taken as na, 7 and the particle 
volume fraction 


(j> = na 3 n/6 , (7) 

where n is the number density of particles. For truly charge-stabilized suspensions, however, the physical hard-core 
is masked by the strong Yukawa repulsion, i.e. the likelihood for two or more spheres being in contact is negligibly 
small. The appropriate physical length scale for these so-called point-Yukawa systems is the geometric mean particle 
distance (r) = n _1 / 3 , and the thermodynamic state and in particular the phase boundaries are determined by two 
parameters only. The phase boundaries of the point-Yukawa system look particularly simple, with nearly straight 
lines, in the two-dimensional (A, T) phase diagram representation [68, 69] where 

A = n(r) , 
f k B T 
u « r » 

are the reduced screening parameter and the inverse reduced Yukawa interaction parameter, respectively. In terms of 
these parameters, the dominating Yukawa-part of the HSY potential reads 

u{x) exp{—A (x — 1)} 
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FIG. 1. Schematic ^A, T^j phase diagram of the hard-core plus repulsive Yukawa (HSY) system with masked hard-core in¬ 
teractions (point-Yukawa system). The arrows indicate the two fluid-state pathways towards the fluid-bcc and fluid-fcc phase 
boundary lines, respectively, followed in our calculations of the short-time rotational self-diffusion coefficient. The left ar¬ 
row corresponds to A = 3 and T = 0.1,0.2,0.5,1,2,5,10,20,50,100,200,500,1000, and the right one to A = 8 and T = 
0.4, 0.6, 0.8,1, 2, 5,10, 20, 50,100, 200, 500,1000. The volume fractions considered in both pathways are (j> = 0.05, 0.15, 0.25, 0.35. 


where x = r/(r) with x > a/(r). If considered as a function of x and q{r), the RDF and static structure factor of 
point-Yukawa particles are uniquely determined by the state point (A, T). 

A sketch of the (A , T) phase diagram of point-Yukawa particles is given in Fig. 1. It consists of a high-temperature 
supercritical fluid phase, separated by a fluid-solid coexistence boundary from a face-centered-cubic (fee) solid phase 
region at high screening, and a body-centered-cubic (bcc) solid phase region at low screening. There is a single triple 
point of three-phase coexistence at A t « 6.9 [ 68 , 70]. According to [71, 72], the fluid-solid coexistence boundary 
determined in simulations is well reproduced by the RY integral equation scheme in conjunction with the Hansen- 
Verlet criterion S(q m ) = 3.1 for the onset of freezing, where q m is the wavenumber position of the static structure 
factor maximum. 

The full phase diagram of the HSY model including systems with significantly non-zero RDF contact values g{r = 
<r + ) > 0 is more complicated, and requires the specification of a third reduced parameter in addition to, say, A and T , 
namely the volume fraction cj). As shown in simulations by Hynninnen and Dijkstra [70], there is then an additional 
triple point at very small A associated with large volume fractions where the fee phase is favored. Provided the 
coupling parameter in Eq. ( 6 ) is sufficiently large (i.e. 7 > 20) and <f> sufficiently small (i.e. <j> < 0.5), the phase 
coexistence lines of the HSY model can be essentially mapped to those of the point-Yukawa model, by expressing 7 
and /ter in Eq. ( 6 ) in terms of {A, X 1 , </>} using (r) oc 0 -1 / 3 . Note that for T > 1, the potential energy of the Yukawa 
potential part at mean particle distance is smaller than the thermal energy ksT . With increasing T and fixed A and 
t f >, the importance of the Yukawa potential part diminishes, and the suspension becomes increasingly hard-sphere like. 

The present study of short-time rotational diffusion is restricted to the fluid phase regime. However, it is interesting 
to compare changes in rotational diffusion when the fluid-bcc and fluid-fcc parts of the fluid-solid coexistence lines 
are approached, respectively, on decreasing the reduced inverse Yukawa interaction parameter T. To this end, in 
our simulation and theoretical calculations of D r we follow two distinct pathways indicated by the two arrows in 
the (A,Tj diagram in Fig. 1. The left pathway is the vertical line along A = 3 with the reduced inverse Yukawa 

interaction parameter series T £ {0.1,0.2, 0.5,1, 2, 5,10, 20, 50,100, 200, 500,1000}, where the smallest value T = 0.1 
describes a state point close to the fluid-bcc phase boundary line part of the point-Yukawa phase diagram. The right 
pathway in the figure is the line along A = 8 with values T £ {0.4,0.6,0.8,1, 2, 5,10,20,50,100, 200, 500,1000}. Here, 
the lowest value T = 0.4 is close to the fluid-fcc phase boundary line part. For both pathways, the volume fraction 
is selected as <f> = 0.05, 0.15, 0.25 and 0.35, respectively. This amounts to simulation-based calculations of D r at 104 
different fluid-phase state points. The static structure factors S(q) and the associated RDFs g(r ) of all HSY systems 
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explored in the present work have been calculated using the RY integral equation scheme described in the following 
subsection. We have checked that each of the considered S(qY s qualifies as a liquid-state structure factor according to 
the empirical Hansen-Verlet criterion. This criterion states that at freezing into a solid phase, S(q m ) attains a value 
near to 3.1 for point-Yukawa systems. For HSY systems with RDF contact values g(u + ) significantly larger than zero, 
the values of S(q m ) at freezing vary in between 3.1 and 2.85, with the lower value attained by a pure hard-sphere 
system [69, 73-75]. 


A. Rogers-Young scheme 


The revised S"f method and the PA scheme require g(r) as the only input. We obtain g[r) numerically by solving 
the Ornstein-Zernike equation [76], 


h{r)=c(r) + n Jd 3 r'c(r')h(\r — r'|) 


( 11 ) 


for a three-dimensional, homogeneous and isotropic fluid in conjunction with the approximate RY [77] closure relation 
invoking the HSY pair potential, 


u(r) 
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( 12 ) 


Here, h(r) = g{r) — 1 is the total correlation function, c(r) is the direct correlation function, and n is the particle number 
density. Eq. (12) includes the mixing function fir) = 1 — exp{—or} with the non-negative inverse length parameter 
a. This parameter is determined self-consistently by requiring equality of the isothermal osmotic compressibilities 
derived from compressibility equation, 
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and the numerically differentiated virial pressure equation, 
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(14) 


where P c and P v is the isothermal osmotic pressure in the compressibility and virial equation, respectively. 

In using Eq. (14), we neglect any thermodynamic state dependence of the pair potential. As noted further up this 
section already, a consequence of integrating out the microionic and solvent degrees of freedom is that the resulting 
effective pair potential of HSY form is in general dependent on the particle concentration n and the system temperature 
T (see, e.g., [66]). This gives rise to additional terms on the right-hand-side of Eq. (14) invoking the partial derivative 
of u(r) with respect to n and T. The precise form of the effective pair potential depends on the specific electro¬ 
chemical surface properties of the colloidal spheres, and the specific properties of the suspending electrolyte solution. 
Since we are not dealing here with the microscopic theory of effective colloidal pair potentials but with the generic 
behavior of D r , we disregard any specific state dependence of u(r), and of the RY mixing parameter a. 

Our numerical solution of the RY integral equation for broad ranges of volume fractions, screening and interaction 
parameters has been facilitated by using a spectral solver described in Ref. [65]. The good accuracy of the RY 
approximation for HSY systems was demonstrated in various studies [21, 40, 52, 53, 78, 79] comprising comparisons 
with simulation and experimental data. 


B. Radial distribution function in the RY scheme 

Owing the rather steep 0(l/r 6 ) long-distance decay of the rotational mobility tensor fj,]]' associated with D r (see 
Sec. VI), the rotational diffusion coefficient is quite sensitive to the shape of the RDF at small particle separations. 

This motivates the following discussion on the behavior of g(r) in the ^A, T, fluid-phase parameter regime of the 
HSY model explored in this work. 

In Fig. 2, the RY calculated RDFs (upper panel) and structure factors (lower panel) for 0 = 0.25 and A = 8 
are depicted for different inverse Yukawa interaction parameters T as indicated in the figure. The strength of the 
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FIG. 2. (Color online) Upper panel: RDF, g(r), predicted by the Rogers-Young scheme for a HSY fluid system with A = 8 
and <j) = 0.25. Various inverse Yukawa interaction parameters T = 0.4,1, 5,1000 are considered as indicated. The horizontal 
dashed line segment marks the Carnahan-Starling (C-S) contact value of hard spheres given by Eq. (15). Lower panel: Static 
structure factor, S ( q ), corresponding to the displayed radial distribution functions in the upper panel. Dashed horizontal line 
segment: C-S compressibility factor of hard spheres according to Eq. (16). Pair distance r and wavenumber q are scaled with 
the geometric mean particle distance (r). 


Yukawa potential at mean particle distance, in units of the thermal energy, decreases with increasing T. For the 
largest considered value T = 1000, the HSY system at <j> = 0.25 reduces essentially to a hard-sphere fluid, with 
the RDF maximum g(r m ) located at contact distance r m = a. With decreasing T, the increasingly strong Yukawa 
repulsion reduces the relative probability, g{cr + ), of two-sphere contact, and it also lowers the compressibility factor 
5(0). Moreover, the nearest neighbor shell of spheres around the radial distance r m where g(r) has its maximum 
moves outwards and sharpens with decreasing T. For the lowest considered value T = 0.4 corresponding to a fluid 
state point near to the fluid-fee phase boundary line, the hard core of the particles is masked and r m ~ (r). 

A measure of the importance of the hard-core part of u(r) relative to the Yukawa part is given by the RDF contact 
value g (<r + ) plotted in Fig. 3 for all considered fluid-phase points (A, T, <(>)■ With increasing T, the relative strength 
of the Yukawa potential part ceases, and a plateau region of the RDF contact value is approached, characteristic for 
hard-sphere-like behavior. This is obviated by the horizontal line segments shown at the right ordinate of the figure 
which mark the Carnahan-Starling (C-S) RDF contact values of hard spheres with vanishing Yukawa tail repulsion 
(7 = 0). given by 


We quote in addition the C-S equation, 


9 C h S s(c 
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FIG. 3. (Color online) Contact value, g(cr + ), of the HSY-RY g(r) as a 
0.05,0.15,0.25,0.35 as indicated. Upper panel: A = 3. Lower panel: A = 8. 
the C-S hard-sphere contact values according to Eq. (15). 


function of T, for volume fractions <j) = 
Horizontal line segments for large T mark 


for the compressibility factor of hard spheres. The origin of the excellent accuracy of the semi-phenomenological 
C-S expressions for hard spheres is still a riddle, and a topic of ongoing research [80]. We emphasize that the 
employed RY scheme is a quite accurate but nevertheless approximate integral equation. Its partial thermodynamic 
self consistency does not imply, e.g., perfect agreement of the RY contact value for hard spheres with the practically 
exact Carnahan-Starling result in Eq. (15). In fact, the RY scheme is lacking thermodynamic self-consistency with 
respect to any thermodynamic property except for the pressure. An extended version of the RY scheme (named 
ERY scheme) has been introduced by Carbajal-Tinoco [ 1], which further improves the accuracy of the original RY 
scheme by introducing a second mixing parameter. In this more elaborate scheme, which however is not applicable 
to pure hard spheres without soft repulsion, the two mixing parameters are determined by enforcing thermodynamic 
self-consistency both regarding the pressure and the excess internal energy per particle. For the sake of simplicity 
and numerical stability, we have refrained from using the ERY scheme in the present work. 

On first sight, it is surprising that in accordance with Fig. 3, the likelihood of observing two closely spaced particles 
is larger for A = 3 than for A = 8, for equal </> and T, even though the screening length of the Yukawa tail is significantly 
shorter in the latter case. This can be understood as follows: While the potential value f3u(x = 1) = 1/T at r = (r) 
is equal for both A values, the repulsive force, —f3du/dx(x = 1) = (1 + A) /T, is larger in the A = 8 case. Taken 
together with the substantially steeper rise of the Yukawa potential with decreasing x for A = 8, this explains the 
lower probability of finding two closely spaced particles. For fixed A, the small-T region where the hard-core is masked 
(i.e. where g(<J + ) ~ 0) shrinks with increasing </>, as it is expected. 

For the upcoming discussion of D r , it is relevant to investigate how the principal RDF maximum g(r m ), and its 
position r m , depend on the pair potential parameters. We notice first from Fig. 2 that r m equals the smallest radial 
distance r > a where the derivative of g{r) turns negative. 

The dependence of g(r m ) and r m on T is shown Fig. 4 for A = 3, and in Fig. 5 for A = 8. According to both figures, 
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FIG. 4. (Color online) Upper panel: Inverse reduced position, {r)/r m , of the principal peak of g(r) as a function of T, for 
A = 3 and values of <j> as indicated. The horizontal line segments at large T indicate the inverse reduced contact distance, 
(r)/a, for the respective cj> values. Lower panel: Principal peak value, g{r m ), of the RDF for the same set of parameters. The 
dashed horizontal line segment at small T indicates the one-component plasma isochoric freezing transition value attained in 
the zero-screening limit A —> 0 (see [72]). 


for fixed <j> and therefore fixed mean particle distance (r), the inverse principal peak location, (r)/r m , increases with 
decreasing Yukawa potential strength, i.e. increasing T, towards the limiting inverse reduced contact distance, ( r)/a , 
of neutral hard spheres. The limiting hard-sphere values for the considered (j> are indicated by the horizontal solid line 
segments at the large-T end of the upper panels in Figs. 4 and 5. Except for A = 8 and the lowest considered volume 
fraction </> = 0.05, the hard-sphere limiting values have been all reached for T = 1000. It is for this (X,4>) point where 
the minimum of g(r m ) as a function of T observed in all depicted curves in the lower panels of Figs. (4) and (5) has 
its largest T value. The minimum in g(r m ) originates from the competition of Yukawa repulsion and excluded volume 
interaction (see here again Fig. 2 for g(r)). This competition is enforced with increasing concentration, reflected by 
a more pronounced minimum moving inwards to smaller values of T. 


We finally notice that the only exception from the monotonic behavior of r m as a function of T is the curve in the 
upper panel of Fig. 4 for A = 3 and (j) = 0.05. To explain the peculiar shape of this curve, in Fig. 6 we plot associated 
RDFs, for values of T as indicated. For large T > 50, the principal maximum of g(r ) is located close to r = er. 
The maximum decreases and is shifted to larger distances r m with decreasing T. At T = 20, no localized principal 
maximum is present any more, and g(r) is monotonically increasing with increasing r so that r m = oo. When T is 
further decreased, the strengthened Yukawa repulsion causes the reappearance of a RDF maximum at a distance r m 
significantly larger than a which is decreasing towards (r). 





















J/<J> 



1 10 100 1000 
T 


FIG. 5. (Color online) Upper and lower panels: Same as in the upper and lower panel of Fig. 4, respectively, but for A = 



bO 



FIG. 6. (Color online) HSY-RY g ( r ) for </> = 0.05 and A = 3. Employed values of T are T = 10, 20, 50,100 as indicated. 
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IV. REVISED BEENAKKER-MAZUR METHOD 

We present here the details of our revised method in its application to rotational self-diffusion. The calculation 
of short-time rotational self-diffusion coefficients on basis of Eq. (5) is not straightforward. The difficulty lies both in 
the calculation of the equilibrium probability distribution function, and of the 3 N x 3 N rotational mobility matrix 
fi rr . The latter difficulty is made explicit when the 3x3 tensor elements of the mobility matrix are represented in 
the form of a hydrodynamic scattering series, 

tij (Ri ■ • • Rjv) = P r Tij (Ri... R n) P r , (17) 

with Tij given by 

Tij (Ri... R w ) = 5ijM (R*) + (1 - 6 io ) M (R,) G (Ri, Rj) M (R,) + 

N 

J2 M (Rj) G (Ri, Rfc) M (Rfc) G (R fc , Rj) M (Rj) + ... . (18) 

k=i, 

The usage of scattering series in the context of suspensions is well established [ L 8 , 82]. Not to interrupt unnecessarily 
our line of reasoning, we therefore refer to the appendix for the explicit expressions for the matrices M and G , and 
the definition of the projection operator P r in Eq. (17). It is worth while to note here, however, that the matrix G is 
related to the Green function tensor, 

G » (r) = 8 ^; (1+ff ‘ ) ' (19) 

of the Stokes equations for an unbounded infinite fluid where r = r/r. This so-called Oseen tensor describes the flow 
generated by a point force at the origin. Moreover, the matrix M (Rj) describes the hydrodynamic response of a 
single sphere with its center at position Rj. Each term in Eq. (18) is called a scattering sequence. Two examples of 
a scattering sequence are 

M (Ri) G (Ri, R 2 ) M (R 2 ) G (R 2 , Ri) M (Ri) G (R 1; R 2 ) M (R 2 ) (20) 

and 

M (R x ) G (Ri, R 3 ) M (R 3 ) G (R 3i R 2 ) M (R 2 ) . ( 21 ) 

Each scattering sequence is a succession of matrices M which scatter the flow, and matrices G which freely propagate 
the flow between scattering events. This physical interpretation of a scattering sequence is useful in our further 
discussion, and for the description of the physical idea underlying the Beenakker-Mazur (BM) method. 


A. Fluctuation expansion 

In one of their first papers on short-time transport properties of suspensions [83] , Beenakker and Mazur introduced 
the so-called fluctuation expansion. The physical idea behind this expansion is related to the resummation of a certain 
class of scattering sequences. BM based their considerations on the scattering series given by the Eq. (2.2) in their 
Ref. [ 13]. In the present analysis, we formulate the fluctuation expansion on basis of the scattering series given by 
Eq. (18). Both approaches are equivalent but the presentation of the BM theory is more straightforward using our 
scattering series representation. 

The scattering sequence in Eq. (21) starts at particle 2 at position R 2 and ends at particle 1 at position Ri, with 
particle 3 acting as the intermediate. The scattering sequence in Eq. (18), with the first propagators G starting from 
the particle 2 and the second one ending at the particle 1, can have all the other N — 2 particles as intermediates, 
summing thus up to the expression 

N 

Y M ( R i) G (Ri, Ri) M (Rj) G (Ri, Ra) M (R 2 ) . (22) 

i—3 

If we treat this expression on the mean-field level by neglecting correlations between the particle positions, it is 
approximated by the volume integral 

n J d 3 R 3 M (Ri) G (Ri, R 3 ) M (R 3 ) G (R 3 , R 2 ) M (R 2 ) 


(23) 
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FIG. 7. Schematic representation of scattering sequences resummed in the fluctuation expansion, a) Scattering sequence in 
Eq. (23). b) Scattering sequence in Eq. (24). c) Scattering sequence with five propagators. 


with the number density n playing in a homogeneous system the role of the single-particle distribution function. The 
procedure of summing up similar scattering sequences involving three propagators G results in 

n 2 J d 3 R 3 j d 3 R A M (Ri) G (Ri,R 3 ) M (R 3 ) G (R 3 , R 4 ) M (R 4 ) G (R 4 , R 2 ) M (R 2 ). (24) 

Notice that in the scattering sequences described by Eqs. (23) and (24), every reflection in the sequence is directed 
towards a new particle. This is schematically shown in Fig. 7. An example of a scattering sequence not satisfying this 
condition is given by Eq. (20). 

To proceed further, it is convenient to follow BM by introducing integral density kernels. Therefore, instead of 
Tij (Ri... Rat) in Eq. (18), we consider its kernel density T defined by 

N 

T(R, R'; Ri.. . Rjv) = ^ 5 (R - R, : ) T {j (R 4 ... Rat) 6 (R' - R,) , (25) 

i,j=1 

where <5(R) is the three-dimensional delta function. The kernel T is represented now as follows, 

T = M + MGM + MGMGM + ... (26) 

where the densities M (R, R'; R 4 ... Rat) and G (R, R') are defined by 

N 

M (R, R'; Ri... R w ) = S (R - R') ^ Af (Ri) <5 (R - R>) , (27) 

i=1 


and 


riR r'i = / 0 for R = R ' 

1 ’ J ( G (R, R') for R/R' ’ 

respectively. In Eq. (26), products of kernels such as MG' appear which should be interpreted as: 


MG 


(R, R'; Ri... Rat) = / dR"M (R, R"; R 4 ... Rat) G (R", R') . 


(28) 


(29) 


Notice that Eq. (26) is a short-hand notation since the kernel variables, the integral signs and the position vectors 
of the particles are omitted. The full notation is used, in contrast, in Eq. (29). It is worth noting that the only 
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difference between G (R, R/) and G (R, R') in Eq. (28) is at R = R'. In introducing G instead of G, BM have 
avoided the summation over the same particles since terms such as M (Rj) G (R,, R, ) M (Rj) are zero. The vanishing 
of G (Rj. Rj) for i = j excludes reflections to the same particle which are absent in the initial scattering series in Eq. 
(18). 

The sequence in Eq. (23) is expressed in terms of the introduced kernels by 


M (Ri) G(M)G (R!,R 2 )M(R 2 ) 


(30) 


and the sequence in Eq. (24) by 


Af(Ri) G(M)G(M)G (Ri,R 2 )M(R 2 ) 


(31) 


As throughout this paper, (...) denotes the average with respect to the equilibrium configurational probability density 
function of N particles. Resummation of all scattering sequences for which each propagator connects to a different 
particle such as in the schematic scattering sequences in Fig. 7, yields 

M (RO G< m) (Ri, Ra) M (R 2 ) . (32) 

Here G^) is an ’’effective propagator”, defined by 

G ( m) =g(i-(M)g)’ 1 , (33) 

which is the sum of a geometric series. The resummation of the above class of scattering sequences in the scattering 
series in Eq. (26) leads to the so-called fluctuation expansion. The resummation can be performed by rewriting 
Eq. (26) as follows 


7 = M + MG 


1 - MG 


M. 


By adding and subtracting (M), and using the operator identity 


(1- A-B)- 1 


{l-A)- 1 


1 — 13(1 — A )" 1 , 


(34) 


(35) 


for two operators A and B , we obtain the fluctuation expansion of the propagator G 


1 - MG 


as 


G 


1 - MG 


-l 


= G 


1- (M- (M))G- (M)G 


-l 


= G 


<M) 


1 — (M— (M» G(m) 


n -l 


When this expression is inserted in Eq. (34), the fluctuation expansion result 


T = M + MG 


(Mi 


1 — (M— (M»G <m> 


n -i 


M, 


(36) 

(37) 


(38) 


for the kernel T is obtained. 


B. Renormalized fluctuation expansion 

The truncation of the fluctuation expansion for the translational mobility matrix has led to an approximate method 
of calculating the short-time translational self-diffusion coefficients [83]. However, the results by this second-order 
fluctuation expansion were found by BM to be unsatisfactory. Therefore, in subsequent work, BM performed another 
resummation which has led to the renormalized fluctuation expansion [13]. 

The scattering sequences illustrated in Fig. 7 for which each propagator links two different particles are resummed 
in the fluctuation expansion. In the renormalized fluctuation expansion, BM resummed similar scattering sequences 
using now a renormalized single-particle operator Mr. Namely, instead of the bare M operators, there appear now 
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FIG. 8. Schematic representation of scattering sequences resummed in a renormalized fluctuation expansion invoking the 
renormalized propagator in Eq. (41). a) Scattering sequence without renormalization of the single-particle response (present 
also in the fluctuation expansion), b) and c) Scattering sequence with the single-particle operator M being renormalized. 


’rings’ built of scattering sequences with renormalized operators. This is illustrated schematically in Fig. 8. BM refer 
to these structures as ’ring self-correlations’, which are scattering structures of the following form 


M fl = M (1 — G 




where the integral operator G? Mh , is defined by 


g \ m r) (R..R0 = 


G <Mr) (R,R') for R = R' 


for R/R' 


while the renormalized propagator G™ s ) is of the form 


g (M r) =g(i — (m*)g)' 


(39) 


(40) 


(41) 


The derivation of the expansion in renormalized density fluctuations is omitted here for it involves only simple 
algebraic manipulations which are given in the works of BM [13-15]. The most important result of this derivation is 
the expansion of the operator 7 in renormalized fluctuations, Mi? — (Mr), according to 


7 — M + MG^h) 

The operator G(m r ) i n the above expansion is defined by 

G 

It is also worth quoting the intermediate result 


1 - (Mr - (M_r,)) G(m k > 


Mi?. 


<m r ) (Ri,R 2 ) — G{m r ) (Ri,R 2 ) — (Ri, R 2 ) • 


(42) 

(43) 


(l - mg) 


G 1 - MG M = G 


<3Vt R > 


1 - (Mi?, - (Mi?)) G (Mr) (Mi?, - (Mi?)) 


G 


<m r ) 


l_(Mi?-(Mi?))G (MR> (Mi?), 


(44) 
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used in the renormalized fluctuation expansion of the translational and rotational self-diffusion coefficients discussed 
in the following. 


C. Renormalized fluctuation expansion for translational self-diffusion 


Before considering the renormalized fluctuation expansion of the rotational self-diffusion coefficient, we consider 
first translational self-diffusion which was investigated earlier by Beenakker and Mazur. The short-time translational 
self-diffusion coefficient D 4 is expressed by the following formula 


D 4 


k B T 

3 


lim 

OO 


Tr 





(45) 


in analogy with Eq. (5) for the rotational self-diffusion coefficient. The above formula has been re-expressed by BM 
in terms of kernels according to 


D t = k ^ p t dpt ^ . 
o 


(46) 


where the matrix d is defined as follows: 


d = nM (R) + M (R) lim 

OO 




M ) (R, R) 


(47) 


The operator P t whose explicit form is given in the appendix projects on the translational components of the matrix 
d. This should be compared with the Eq. (3.16) in [83]. The operator expression in Eq. (44) inserted into the above 
expression results in the renormalized fluctuation expansion of d. Up to second order in the renormalized fluctuations, 
Mr — (M a) , this expansion reads 


where 


and 


d = d<°> + dP> + d P) + ..., 


d (0) = nM + M [G (Mfl> (Mr)] (R, R) , 
d (1) = 0, 

d^ = d^ + d g 2 ), 



M 

M 


G 


(Mr) 


G 


(Mr) 


((Mr - (Mr))G {M r} (Mr 
{(Mr - (M R ))G {Mr} (Mr 


W) 

(Mfl)) 



) ^<M r ) (Mfl) 


(R,R) . 


(48) 


(49) 

(50) 

(51) 

(52) 

(53) 


The truncation approximation, 


d « d (0) + d (1) + d (2) , (54) 

constitutes along with Eq. (46) and Eqs. (48-53) the second-order renormalized fluctuation expansion approximation 
for the short-time translational self-diffusion coefficient. 


D. Renormalized fluctuation expansion for rotational self-diffusion 

The extension of the renormalized fluctuation expansion method (5"f scheme) to short-time rotational self-diffusion 
was made first by Treloar and Masters [37]. This extension is straightforward in our formalism, since Eq. (5) for 
D r is similar in structure as Eq. (46) for D 4 . The only difference appears in the invoked projectors, i.e. instead of 
the projector P 4 on the translational components the projector P r on the rotational components is used in Eq. (5). 
Since the renormalized fluctuation expansion has been introduced in the previous section on the general level of the 
d matrix, the extension from translational to rotational self-diffusion is straightforward. 
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E. High-accuracy second-order renormalized fluctuation expansion 

Next, we point out the differences between the original second-order BM method of calculating D r by Treloar 
and Masters, and our revised second-order dy scheme. Approximations in the dy scheme are made in particular 
in two calculation steps. The first one is the truncation of the series in renormalized fluctuations to second order 
described in Eq. (54). Secondly, also the matrices S 2 ' 1 are approximated since they relate to infinite dimensional 
hydrodynamic matrices such as M, G, Gy yt H >, or G s , Mr ^ which for a numerical evaluation must be truncated (see the 
appendix). In the works by BM on the hydrodynamic function and the high-frequency effective viscosity, and in the 
work by Treloar and Masters on rotational diffusion, due to technical difficulties severe truncations of these matrices 
have been made as discussed in detail in Ref. [24]. In our revised method, we also truncate the infinite dimensional 
hydrodynamic matrices, but different from earlier works an extrapolation to infinite dimension has been included. 
The details of this extrapolation procedure are along the lines described in Ref. [2 ] and are thus not repeated here. 


V. SIMULATION METHOD 

We have calculated D r to high precision for no-slip spheres using a hydrodynamic multipole method corrected for 
lubrication [36, 84-86], and encoded in the HYDROMULTIPOLE program package [36]. The values for D r have 
been determined from equilibrium configuration averages using typically N = 256 spheres interacting with the HSY 
potential, and placed in a periodically replicated cubic simulation box. At least 150 independent configurations for 

each parameter set ^A, T, <pj were used. This has resulted in a statistical relative error of less than 0.001. As reported 

in Ref. [32], the calculated values for D r (N) using the periodic simulation box with N particles are not critically 
dependent on the system size. Therefore, no system size correction extrapolating to an infinitely large system is 
required as for short-time collective diffusion properties. 


VI. PAIRWISE ADDITIVITY APPROXIMATION 


The rotational hydrodynamic mobility tensor //V(Ri... Rat) of N spherical particles in an infinite, quiescent fluid 
linearly relates the hydrodynamic torque T, acting on a particle i to its rotational velocity By disregarding 
three-body and higher order hydrodynamic cluster contributions, on can approximate the exact iV-particle rotational 
hydrodynamic mobility tensor by a sum of two-particle contributions, 


/xjaRi-.-Rjv) 


r 

Mo 


N 

i+ y , w u( R ‘ 


R-n) 


(55) 


where /Iq = D^/ksT is the single-sphere rotational mobility coefficient and 1 the three-dimensional unit tensor. The 
two-sphere tensor (R, — R n ) describes the hydrodynamic self-interaction of sphere i by means of flow reflections at 
a second sphere labeled by n, in the absence of the N — 2 other particles. This constitutes the pairwise additivity (PA) 
approximation where it is assumed that the His between two spheres are not disturbed by other ones. In principle, 
this assumption is justified for semi-dilute systems only. 

On exploiting the axial symmetry of the two-sphere problem, the two-sphere tensor can be split into longitudinal 
and transversal components, 


«E(Ri - Rn) = aK(Rin)&inRin + P[[(R in ) 


1 R,„R 




(56) 


with R.y 77 — (R-y R’n) / A -,77 and R%n — Ry R/i | ■ 

In term of the scalar longitudinal and transversal functions aYi{R) and /3ii(R), the normalized short-time rotational 
self-diffusion coefficient is expressed in PA approximation by [3, 38, 39] 


D r 

Rn 


oo 

= 1 + 8^ J dx x 2 g(x) \a r Yi{x) + 2f$[\{x)] 


(57) 


where x = r/a is the two-sphere distance in units of the sphere diameter cr = 2a. The functions a^(r) and fd{[(r) 
can be calculated recursively in the form of a power series in the reduced inverse pair distance a/r. For the no-slip 
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<|)=0.05 



T 


FIG. 9. Normalized rotational self-diffusion coefficient, D r /D q, as a function of T, for A = 3 and 8, and </> = 0.05. Simulation 
data (plus symbols and crosses) are compared with revised <Sy method predictions (solid and short-dashed curves) and PA 
method predictions (dashed-dotted and long-dashed curves). The theory and simulation values of D r for A = 8 are in general 
larger than the respective ones for A = 3. Horizontal solid segment at large T: Hard-sphere value according to Eq. (1). 
Horizontal solid line segment at small T: Scaling prediction in Eq. (2) for low-salinity charge-stabilized systems. 


hydrodynamic surface boundary condition employed in this work, the leading-order (far-held) contributions are 



(58) 

(59) 


with higher-order terms in the expansion given in [44, 87]. At near-contact distance r ss 2a where lubrication comes 
into play, the expansions in Eqs. (58) and (59) converge only slowly. In our numerical implementation of the PA 
method, we therefore use high-order series expansion results obtained by Jeffrey and Onishi [ ]. On using the 
zero concentration lrard-sphere RDF, g(x) = Q(x — 1), in Eq. (57), where 0(a;) is the unit step function, we have 
numerically checked that our code precisely reproduces the first-order virial coefficient value 0.631 in Eq. (1). This 
demonstrates the high accuracy of the employed tabulated values for a\\(r) and also at near-contact distances. 


VII. RESULTS AND DISCUSSION 

Our simulation and theoretical results for D r /D q as a function of the inverse reduced Yukawa interaction parameter 
T are depicted in Figs. 9 - 12, for volume fractions </> = 0.05 — 0.35. In each figure, the results for the vertical fluid- 
phase pathway at A = 3 directed towards the fluid-bcc phase coexistence line is compared with the results for the 
pathway at A = 8 directed towards the fluid-fcc coexistence line (cf. Fig. 1). Note the different ordinate scales in the 
four figures selected to highlight the differences in the theoretical and simulation results for D r . 

We start with the discussion of the simulation results for the HSY systems. The slowing influence of the His 
on rotational self-diffusion increases when the likelihood for near-distant particle pairs increases. According to our 
discussion of the RDFs in Fig. 3, the particles for A = 8 repel each other more strongly than those for A = 3, so 
that the radial region where g(r ) is small is more extended in the former case. This explains why all curves of D r for 
A = 8 are located above those for A = 3, for all considered volume fractions. The HSY particles can approach each 
other more closely with increasing T. This is reflected in curves for D r which are monotonically decreasing. 

As discussed earlier, at large values of T a plateau region of D r is reached where the particles behave essentially 
as neutral hard spheres, independent of A. The simulation curves in Figs. 9-12 converge therefore for large T towards 
the result in Eq. (1) which accurately describes the ^-dependence of D r for neutral hard spheres up to the freezing 
transition volume fraction. With increasing (j>, the hard-sphere plateau region is reached for smaller values of T. In 
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FIG. 10. Same as in Fig. 9 but for <f) = 0.15. 


the opposite limit of low T values, the interaction of the HSY particles is dominated by the Yukawa potential part. 
For smaller volume fractions where g(cr + ) « 0 is observed, Eq. (2) derived originally for low-salinity charge-stabilized 
particles is expected to be a decent description of D r in the small-T region. It is noticed that the simulation curves 
for A = 3 and 8 , and (j) = 0.05 and 0.15, are indeed converging, with decreasing T, towards the result in Eq. (2). 
Differences are visible for <j> = 0.25 and 0.35 where Eq. (2) provides only an upper bound for D r . 

The simulation curves for both values of A are well reproduced by the PA approximation for <f> = 0.05 and 0.15. 
This is an expected feature of this method which becomes exact at low concentrations. At higher concentrations, 
three-body and higher-order His contributions come into play which are disregarded in the simple PA treatment. As 
a consequence, the simulation curves for <f> = 0.25 and 0.35 are underestimated at all values of T, i.e. the slowing 
influence of the His on rotational self-diffusion is overestimated. As it was noticed earlier in the context of translational 
self-diffusion [51], this can be attributed to the fact that the PA approximation neglects shielding of the His between 
two particles by other particles in their vicinity. While D r at larger <f> is underestimated in the PA approximation, 
the T region where the simulation curves for A = 3 and 8 converge is still well predicted. The PA method for D r 
is actually a decent approximation up to surprisingly large volume fractions. This should be contrasted with its 
performance for collective diffusion properties where its concentration range of application is significantly smaller [9]. 

We discuss now the performance of the revised second-order S'y scheme. It is noticed from Figs. 9-12 that it 
systematically, and significantly, underestimates D r , for all values of T and all considered volume fractions. The 
relative error, | D £ 7 — -D£ im |/-Dg im , increases systematically with increasing <j>, and it is more pronounced at the lower- 
T side where the Yukawa repulsion is strong. In the hard-sphere-like interaction regime of large T values, the relative 
error increases from about 1% at <fi = 0.05 to 23% at <p = 0.35. In comparison, the relative error in the small T region 
is larger, increasing from about 2% at <j> = 0.05 to 32% at <f> = 0.35. The most significant feature of the (revised 
and non-revised) S'y method at higher concentrations is its weak sensitivity to changes in range and strength of the 
pair potential, and to the accompanying changes in the RDF. For example, in Fig. 12 with <j> = 0.35, the relative 
difference between the curves for A = 3 and A = 8 is four times larger for the simulation data than for the S'y method 
curves. Incidentally, a similarly weak dependence on the shape of the pair potential and RDF has been found for the 
non-revised S'y method result for the high-frequency viscosity of charge-stabilized suspensions [9]. 

The revised S'y scheme for D r performs somewhat better for neutral hard spheres. A direct comparison with 
simulation results for hard spheres (where A = oo or T = oo) is made in Fig. 13. The simulation data are well 
described by Eq. (1) in the full liquid-phase concentration range . While the revised S'y method for D r significantly 
improves the original second-order S'y method results for hard spheres by Treloar and Masters in the range < 0.4, 
there is no improvement at larger volume fractions. A general observation made for the HSY systems is that the 
relative mean difference between revised and non-revised second-order S'y results for D r is typically 5 — 7% or less. 

The low sensitivity of the S'y scheme on the shape of the RDF can be related to its mean-field structure. An 
important ingredient of BM method is the effective propagator Gp y[ R ), defined by Eq. (41), which depends on $ 
but not on the RDF or higher-order static correlation functions. The suspension microstructure for a given volume 
fraction enters into the BM approach only through fluctuations which are included up to second order. The truncation 
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FIG. 11. Same as in Fig. 9 but for (f> = 0.25. 


4>=0.35 



T 


FIG. 12. Same as in Fig. 9 but for (f) = 0.35. 


of the renormalized fluctuation expansion at higher order than the second one could arguably enlarge the sensitivity 
of the method on the equilibrium suspension microstructure, for the price that triplet or even higher-order static 
distribution functions are then required as additional input. 

Fig. 13 includes also the PA approximation prediction for hard spheres. The simulation data are well described by 
this method for 4> < 0.2, but D r is increasingly underestimated at larger </>. 


VIII. SUMMARY AND CONCLUSIONS 

We have presented the first comprehensive theory-simulation study of short-time rotational diffusion in suspensions 
in the fluid-like phase, with the particles interacting by the HSY potential. Since this effective pair potential is generic 
to many different soft matter systems including ionic microgels and globular protein solutions, the presented results 
should be of broad interest. 

A large body of high-precision simulation data was generated and compared with the results by two theoretical 
methods. The first and more elaborate one is a revised second-order version of the original Beenakker-Mazur method 
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FIG. 13. Normalized rotational self-diffusion coefficient D r jDg of neutral no-slip hard spheres as a function of <j>. Simulation 
results by Abade et al. [32] and Banchio et al. [21] are compared to the second-order virial expansion result in Eq. (1) by 
Cichocki et al. [36], the original A 7 method result by Treloar and Masters [37], and our revised second-order Ay method and 
PA predictions. 


which has been adapted to rotational diffusion by Treloar and Masters [37], in the context of no-slip hard spheres. 
In our revised second-order Ay method, various approximation steps made in the original method have been avoided, 
in particular regarding the treatment of the His which in the original method was rather approximate. The second 
method is the PA approximation with full account of two-body His contributions including lubrication terms, but 
with three-body and higher-order His contributions disregarded. 

General features of D r observed in our simulation study are its monotonic decrease with increasing T, reproduced 
qualitatively by both theoretical methods, and its strong sensitivity on the Yukawa potential range parameter A for 
intermediate values of T. This sensitivity is well captured by the PA method, different from the revised Ay method 
which shows this sensitivity for low concentrations only. A lower bound of D r is provided by the Eq. (1) for hard 
spheres, reached by the simulation curves of D r at large values of T. An upper bound is given by the scaling result in 
Eq. (2). This bound is approached by the simulation curves at low values of T, provided </> is sufficiently small (i.e., 
(j) < 0.15) and A not very large (i.e., A < 8 ). 

Even though His are accounted for to significantly higher accuracy than in BM and Treloar and Masters approach, 
the resulting improvement of D r in our revised second-order Ay method is comparatively small, amounting roughly 
to 5 — 7% for (j) < 0.4. A similar observation has been made regarding the hydrodynamic function, the short-time 
translational self-diffusion coefficient, and the high-frequency viscosity of hard spheres [24]. This may be due to the 
interplay of mean-field and His approximations going into the A 7 method which can cause uncontrolled fortuitous 
cancellations or fortifications of errors. 

The (revised) second-order £7 method performs distinctly better for collective than self-diffusion properties, namely 
for the wavenumber dependent distinct part of the hydrodynamic function and the collective diffusion coefficient [9]. 
In particular regarding the distinct hydrodynamic function part, the Ay method performs quite well both for hard 
spheres and charge-stabilized particles. In future work, it will be interesting to find out whether the overall good 
agreement with simulation data for the collective diffusion properties of HSY systems can be further improved by the 
revised method. The performance of the revised £7 method regarding D r can be possibly improved by the inclusion 
of third-order renormalized fluctuation contributions where, however, static three-body distribution functions are 
required as input to the extended method in addition to g(r). 

The PA method with its full account of two-body His contributions describes the HSY simulation data for D r quite 
well for volume fractions up to ^ ~ 0 . 2 . At larger concentrations, however, the rotational self-diffusion coefficient 
is underestimated owing to the disregarded hydrodynamic shielding effect embodied in three-body and higher-order 
hydrodynamic mobility tensor contributions. In going beyond the PA approximation, three-body irreducible hydro- 
dynamic cluster contributions could be additionally considered. For low-salinity systems, the leading-order far-distant 
three-body contributions have been accounted for in [38], in conjunction with Kirkwood’s superposition approximation 
for the static three-body distribution function. The long-distance three-body cluster contribution to D r is positive 
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valued, with the effect of bringing the value for D r thus closer to the simulation data [27, 38]. 

Finally, we note that both the revised method and the PA approach can be rather straightforwardly extended 
to colloidal particles with internal hydrodynamic structure, and hydrodynamic surface boundary conditions different 
from the no-slip one used in the present work. This offers the possibility to study theoretically, e.g., the rotational 
self-diffusion of weakly crosslinked ionic and non-ionic microgels, and of core-shell particles with a fluid-permeable 
soft shell. 
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APPENDIX: HYDRODYNAMIC MATRICES 


The translational and rotational mobility matrices with tensor elements /x// (Ri ... Rjv) and /x/( (Ri ... Rjv) ap¬ 
pearing in Eqs. (5) and (45) for the short-time rotational and translational self-diffusion coefficients, respectively, can 
be calculated using the hydrodynamic multipole matrices Z 0 , ho, Z 0 , and G introduced in Ref. [86, 89]. Each of these 
matrices is indexed by the set of three indices (Z, to, a) with 1 = 1,..., oo, m = — l, —l + l,...,l, and er = 0,1, 2. The 
matrix Zq is given by the formula 




“21(TO’ —1 


Zl,oo' 


(60) 


where the dimensionless coefficients zi }ITC ,' have been defined in Ref. [ ]. The only non-vanishing elements of the 

matrix /xq are for l = V = 1 and a = a' = 0, with 


[/-Xo(Ri)]zmO,/ / m / 0 = 6 U >6 77 


9 rja ’ 


and for l = l' = 1 and a = a' = 1, with 


[MO — $ 11 ' 


6rja 3 


(61) 


(62) 


The matrix Zq is related to /xo and Zq by 

^o(R-i) = Zo(Ri) — Zq (Rx)/xo(Rx) Zq (Rj). (63) 

The Oseen tensor in multipole space, G(R,, R ; ), is for non-overlapping configurations given by 

^(RijRjjJzmCT j/mv/ = Ulm 5 ,+ ~(Ri - R j,lmcr, I'm 1 cr 1 ) for |R; - Rj| >2 a, (64) 

Wl'm ' 

where the coefficients S A and ni m have been introduced in Ref. [90] . The matrices noted above are used to construct 
the scattering series for the generalized multipole mobility matrix elements according to 

Hij (R-i • • ■ I Uv) = fiij/J-o (R-i) + (1 — dij) no (Rj) Z 0 (Ri) G (Ri, Rj) Zq (Rj) ho (R-j) + 

N 

+ ^2 Mo (Ri) Zq (Ri) G (Ri, Rjt) Zq (Rfe) G (Rfe, Rj) Z 0 (Rj) ho (Rj) 

k= 1 , 


(65) 
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In the multipole basis, the lowest multipole elements l = l,er = 0 correspond to translational motion. Therefore, 
the tensorial elements of the translational mobility matrix can be calculated from the generalized mobility matrix 

with elements by an appropriate projection. To this end, the projector P is introduced by 


[P*] 


,mcr 



( 66 ) 


(R) 

where a = 1, 2, 3 denote three Cartesian components. The tensor yi m = Xy\ J is defined in Ref. [91], together with 

X and y( R \ Moreover, P* is the Hermitean conjugate of p\ Finally, the Cartesian mobility matrix is related to the 
generalized multipole mobility matrix by the projection operation 


/i« (Ri • ■ • R-Jv) = P'lHj (Ri • ■ • Rat) ■ (67) 

The according expression for the rotational mobility matrix reads 

/4J (Ri ... Rat) = P r ', Mj (Ri... Rjv) R . (68) 

Here, the operator P projects from the multipole space the elements l = l,er = 1 to the Cartesian components 
corresponding to rotational motion, i.e. 


[P r 


\at,lmcr 



(69) 


In the scattering series for the mobility matrix in Eq. (65), the single particle matrices /j.qZq and Zq^lq at the start 
and end of a scattering sequence are different from the single-particle matrices Zq dispersed in between the propagators 
G. In the derivation of the renormalized fluctuation expansion, it is thus very useful to rewrite the scattering series 
in a way that all the single particles matrices starting a scattering sequence or appearing in between the propagators 
G are the same. 

This can be achieved by introducing the matrix defined by [92] 


M = 


Mo Mo-^o 
Zo^O Zq 


(70) 


According to the above definition, the matrix M has in addition to the multipole indexes l,m,a the index u = 1,2. 
Therefore, = [poZo\ima,vm'a' et cetera. Similarly, we define the matrix 

\P*\uln-ia ,u' l'm' a ' by 


0 0 
0 G ’ 


_ f _ l ■ 

and generalize the projectors P and P to 


\ot,ulmcr ^ul^ll^aO \ A [yim]o: > 

m 47j- 


(71) 


(72) 


[P ]a,ulmc t — 



(73) 


In the extended space expressions in Eqs. (65) and (68), the rotational mobility matrix can be rewritten in the form 
given by Eq. (17). This can be done analogously for the translational mobility matrix, with the only difference lying 
in the different projection operators. 
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